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GRAVITATIONAL FRAGMENTATION OF EXPANDING SHELLS. I. LINEAR ANALYSIS 

Kazunari iwasaki^ Shu-ichiro Inutsuka\ and Toru Tsuribe^ 
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ABSTRACT 

We perform a linear perturbation analysis of expanding shells driven by expansions of H// regions. 
The ambient gas is assumed to be uniform. As an unperturbed state, we develop a semi-analytic 
method for deriving the time evolution of the density profile across the thickness. It is found that 
the time evolution of the density profile can be divided into three evolutionary phases, deceleration- 
dominated, intermediate, and self-gravity-dominated phases. The density peak moves relatively from 
the shock front to the contact discontinuity as the shell expands. We perform a linear analysis taking 
into account the asymmetric density profile obtained by the semi-analytic method, and imposing the 
boundary conditions for the shock front and the contact discontinuity while the evolutionary effect of 
the shell is neglected. It is found that the growth rate is enhanced compared with the previous studies 
based on the thin-shell approximation. This is due to the boundary effect of the contact discontinuity 
and asymmetric density profile that were not taken into account in previous works. 
Subject headings: HII regions - hydrodynamics - instabilities - shock waves - stars: formation 



1. INTRODUCTION 

Expanding shells are ubiquitous in the interstellar 
medium. They are driven by energetic phenomena of 
massive stars, such as emission of ionizing photons, 
stellar winds, and supernova explosio ns. Recently, us- 
ing 1 02 samples identified as shell, iDeharveng et al.~l 
(|201Cl l found evidences of the star formation in more 
than a quarter of the shells, suggesting that the trig- 
gered star formation by HII regions may be an effi- 
cient process of the mass i ve sta r formation. Theoreti- 
cally, lElmegreen fc Ladal ()1977|) presented a sequential 
star formation scenario where the massive star forma- 
tion takes place through gravitational fragmentation of 
the expanding shell that is driven by HII regions sur- 
rounding massive stars, and newly formed massive star 
also triggers the formation of next generation. 

To understand the triggered star formation, it 
is important to investigate how and when the ex- 
panding shell fragments through the gravitational in- 
stability (GI). Earliest studies were done by using 
linear analyses of the static dense gas layer con- 
fined by the same th ermal pressures of hot rarefied 
gases on both sides (iGoldreich fc Lvnden-Bell 119651 : 



lElmegreen fc Elmegreen 1119781: iLubow fc PringielllQQSD . 

They showed that the GI begins to develop with a scale 
comparable to the layer thickness and with a growing 
timescale comparable to the free-fall time of the layer. 
However, their linear analyses are oversimplified because 
the actual shells are confined by the shock front (SF) on 
the leading surface and the contact discontinuity (CD), 
or the ionization front (IF) on the trailing surface. More- 
over, an unbalance between the ram pressure and the 
thermal pressure causes a decelerating or an accelerat- 
ing expansion. Many authors have tackled the stability 
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analyses with these effects by mainly using the thin-shell 
approximation where the perturbed variables are aver- 
aged across the thickness. The stabili ty analysis of ex - 
panding shells has been investigated by Vishniac (19831) , 
Elmegreen (1994), and Whitworth et al. (1994b). They 
took into account dilution effects of perturbations ow- 
ing to the expansion and the mass accretion. Their lin- 
ear analyses of expanding shells neglected the structure 
across the thickness and the boundary effect of the CD. 
Thus, how these effects that t hey ne g lected infiuence the 
GI have been unknown yet. iVoit I (|1988) investigated 
stability of asymmetric layers, and found that the asym- 
metry of the density profile of the shell greatly influences 
the development of the GI. Moreover, by using shock-like 
boundary conditions, he found that the different choice 
of the boundary condition greatly modifies the dispersion 
relation. However, their analysis is limited to be in the 
incompressible fluid. 

In this paper, we perform a linear analysis taking into 
account the structure across the thickness and the effects 
of boundaries, i.e., the SF on the leading surface and the 
CD on the trailing surface. In order to determine the 
density profile all the time, we develop a semi-analytic 
method that well describes the one-dime nsiona l (ID) 
evolution. This paper extends the study of iVoit I 



to include the compressible effect and the more realis- 
tic density profile by taking into account the ra dial self- 
gravitational force (jWhitworth fc Francis |[200^ . We ne- 
glect the effects of expansion and mass accretion through 
the SF. 

In this paper, since we focus on investigation of how 
the boundary effects and asymmetric density profiles in- 
fluence the GI, we do not apply our result to estimate 
fragmentation time and scale. We will perform three- 
dimensional simulation of expanding shells to compare 
with the results of the linear analysis, and present de- 
tailed quantitative aspects of the fragme ntation process 
of exp anding shells in a subsequent paper (jlwasaki et al. I 
1201 IL submitted). 

The outline of the paper is as follows: in Section [2j we 
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present a thin-shell model of the expanding shell driven 
by the HII region. In Section [3l we develop a semi- 
analytic method to derive time evolution of the density 
profile. We investigate influences of the asymmetric den- 
sity profile on the dispersion relation of the GI by con- 
sidering pressure-confined layer in Section |4l In Section 
[5l we perform linear analysis of expanding shells by us- 
ing density profile obtained in Section [3] and by imposing 
the approximate SF and the CD boundary conditions. In 
Section |6l we compare our results with previous works. 
Summary is presented in Section [T] 

2. THIN-SHELL MODEL DRIVEN BY HII REGION 

Massive stars emit ultraviolet photons {hv > 13.6 eV) 
and produce HII regions around them. Here, we con- 
sider a massive star that emits ionizing photons with the 
photon number luminosity Quv [s~^], into the ambient 
gas with the uniform density of pE = mnE, where tie 
and m are the number density and the mean mass of the 
ambient gas particle, res pectively. In the standard pic- 
ture (e.g., Spitzer Ill978i ). the IF initially expands with 
a supersonic speed with respect to the sound speed of 
ionized gas, cn. The HII region begins to expand by the 
pressure difference between the HII region and the ambi- 
ent gas when the IF reaches the Stromgren radius, i?sT 
given by 



RsT — 



3Quv 



1/3 



(1) 



where indicate the case-B recombination coefficient. 
In this phase, the SF emerges in front of the IF and 
sweeps up the ambient gas into a dense shell. This pa- 
per focuses on the evolution of the shell after the shock 
emerges. The equation of motion of the shell is given by 



dt r^^ dt 



(2) 



where Ms — AttGpeR^/?> is the total mass of the shell, 
i.e. the mass of the ambient gas that initially occupied 
the volume of the HII region, Rg is the mean radius of the 
shell and Pn is the thermal pressure of the HII region. 
Here, we neglect the pressure of the ambient gas and the 
thickness of the shell. In the HII region, the detailed 
balance between the recombination and the ionization is 
approximately established all the time. Therefore, Pn 
can be expressed using _Rs as follows: 



Pn = PEC„f— j 



3/2 



(3) 



Using Equation ([3]) , we obtain the solution of Equation 



Rsit)^RsT 1 



7 



cut \ 



4/7 



(4) 



/URstJ 

(jHosokawa fc Inutsuka1l2006f l. Equations Q and (g]) are 
valid only in the early phase. As the shell sweeps up the 
ambient gas and increases its mass, the self-gravity infiu- 
ences the expansion. The equation of motion including 
the self-gravity becomes 



dt\ dt 



ATiRiPii 



2i?2 



(5) 



where the second term on the right-hand side represent s 
the self-gravitational force (iWhitworth fc Francis |[200^ . 
The factor of 1/2 in the self-gravity term arises because 
the gravitational acceleration vanishes at the inner sur- 
face, it is GMs/Rg at the outer surface, and the mass- 
weighted average across the thickness is GMs/2R'^. One 
can see that the self-gravity slows the expansion in Equa- 
tion ([5]). 

In this paper, for convenience, the units of the time, 
length, and mass scales are taken to be 



in = 



32Gpi 



= 1.6 n 



-1/2 
E,3 



Myr, 



Rf) — 



'7 + \ 4/7 
ICuto \ „3/7 



R7rr = 5.9 Q 



1/7 



UV,49 "E.a 



4/7 

c pc. 



and 



Mo = peRI = 5.0 X 10^ Q^/^ 



Tip,, AIq, 



UV,49 "E,3 
49 „-l 



(6) 
(7) 

(8) 



respectively, where Quv,49 = Quv/10 s , and tie, 3 — 
tie/IO^ cm~^. 

Non-dimensional quantities normalized by tg, Rq, and 
Alf) are expressed by using tilde, e.g., Rs = Rs/Ro- Using 
non-dimensional quantities, we can rewrite Equations ([3]) 
and (HI) as 



and 



dt I ' dt 




(9) 
(10) 



respectively. We integrate Equation pUj) with respect to 




Fig. 1. — Expansion laws of shells. The abscissa and ordinate 
axes indicate the time t/to and the radius of the shell Rb/Rq, 
respectively. The solid lines correspond to the case with (riE/cm-^, 
Quv/s-i)= (io3^ lo-iS), (102, 10«), (10*, 10«), (lO^, lO**), and 
(103, 1045). 

time with the initial condition, R — _Rst at t = 0. The 
initial velocity di?sT/dt is obtained from Equation (jlj 
with t — 0. Figure [1] shows the obtained expansion law 
with various parameters, (nE/cm~a^ Q\jv/s~^) = (10^, 
10*9), (102, lO^^), (10^ 10^9), (10^, lO^S), and (10^, 
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2q45^_ The difference of these parameter gives the dif- 
ferent values oi Rs at t = as shown in Figure [TJ In 
Figure \T\ it is seen that as the shell expands, the lines 
quickly approach to an asymptotic line that is indepen- 
dent of the parameters. Therefore, the dependence of 
the expansion law on the parameters is approximately 
el iminated by using the non-di mensional quantities. 

iWhitworth fc Francis I ()2002D derived similar thin-shell 
equations for the shells driven by steady stellar winds. 
They found a change of the power- law index (from 3/5 to 
1/5) at the time when the self-gravity starts to be impor- 
tant. In the gravity dominated phase, the shell expands 
keeping the force balance between the thermal pressure 
of the hot bubble and the self-gravity. In the stellar wind 
case, the steady energy input allows the outward expan- 
sion of the shell (cx rI^^) even when the self-gravity be- 
comes important. On the other hand, in the case with the 
HII regions, from Equation ^TU\\ . the gravitational force 
(cx i?g ) increases more rapidly than the pressure force by 

"1/2 

the HII region (cx i?s ), suggesting that the shell begins 
to collapse toward the center at a certain radius. In the 
numerical calculation, it occurs at i?s ^ 2.3 in all param- 
eters. The last term of Equation ([TOl) is valid until only 
the expansion phase. In reality, besides ionizing photon, 
the massive star emits strong stellar wind continuously 
over several tens of million years and dies through su- 
pernova explosion (Weaver et alTI 119771 ). They may in- 
fluence the dynamics of the shell in the self-gravity dom- 
inated phase. In this paper, for simplicity, we focus on 
the expansion phase by the ionizing photon. 

3. TIME EVOLUTION OF DENSITY PROFILES: 
UNPERTURBED STATE 

In this section, we derive the time evolution of the den- 
sity profile of the shell in a semi-analytic way. We assume 
that the shell is in instantaneous hydrostatic equilibrium 
at each instant of time. This is reasonable assumption 
because the shell is very thin and the sound-crossing time 
across the thickness is very short compared with the ex- 
pansion timescale. The equation of the hydrostatic equi- 
librium in the frame of the shell is given by 



Cgdp d0 
p ar dr 



0, 



(11) 



where gdcc = —d^Rs/dt^ is the inertial force owing to 
the deceleration of the shell, and is assumed to be spa- 
tially constant within the shell. In the decelerating shell, 
the inertia force is parallel to the radial direction. The 
Poisson equation is 



d20 2d0 d^c/. 

TT + ~:r - TT = ^^Gp, 
dr r dr dr^ 



(12) 



where the curvature effect is neglected because Rs is 
much larger than the thickness. We confirmed that the 
curvature effect is negligible by comparing density pro- 
files with and without curvature effect. Substituting 
Equation (|lip into Equation ([T^ . one obtains 



dr ( p dr j ^ 



(13) 



Equation (|13l) can be solved analytically as follows: 

Pir) = POO jcosh (^^^) } > (14) 

where Rc and poo are the ra dius and the d ensity where 
dp/dr — 0, respectively (c.f. ISpitzer lll942[ ). and Hq = 
Cs/V2ttGpoo is the scale height. 

From Equation ()14[) . if we determine poo and Rc, 
the density profile is completely specified except for the 
boundaries that are discussed later. Here, the value of 
Rc itself loses its physical meaning since the curvature is 
neglected. Therefore, only pqq specifies the density pro- 
file. The peak density poo is determined by the condition 
of the force balance at r = i?cD • The gravitational force 
must vanish at r = i?CD because the total mass of the 
hot bubble is negligible. Therefore, from Equation dTTT) . 
the inner boundary conditions are given by 



c|dp 
p dr 



5do 



(15) 



The column density from i?cD to Rc is 



Sdoc — 



pdr = pooHq tanh 



Rc — RcB \ 



.9do 



47rG" 
(16) 

where we use Equation (jl5p in the last equality. The 
characteristic column density Sdoc represents the amount 
of the deceleration. The ratio of the column density Sg 
to Sdcc determines the importance of self-gravity relative 
to deceleration. From Equation and the pressure 
equilibrium at the CD, p(i?cD)Cg = Pu, the peak density 
can be expressed by Sdcc and Pu as follows: 

Pooc^ - ^11 + 2^GSL- (17) 

Substituting Equation ([T5| into Equation ([T7|) . one ob- 
tains 



5dc 



SttG 



(18) 



The peak density poo is determined by the following way. 
We use the thin-shell model shown in Section [2] to get 
d^Rs/dt^ = —gdcc and Rs at any given times. The pres- 
sure of the HII region Pu is given by Equation ^ . Sub- 
stituting obtained gdcc and Pu into Equation (jlSp . we 
can get poo, and can specify the functional form of the 
density profile. 

Next, we determine the positions of boundaries, both 
the CD and the SF. As mentioned above, only the dis- 
tance relative to Rc has physical meaning. The position 
of the CD and the SF is determined from the pressure 
balances on both sides which are given by 



c^piRcD)^Pn and CsP{Rsf) = Pe 



dRs 



(19) 



respectively, where dRs/dt is obtained from the thin-shell 
model. 

3.1. Three Evolutionary Phases of Density Profiles 

The time evolution of the density profile is character- 
ized by the ratio of the column density Sg to the char- 
acteristic column density Sdec- Here, the column den- 
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a) deceteration-dominated phase 




(b) intermediate pliase 



PCD 




c) self-gravity-dominated pliase 



PCD 




Fig. 2. — Schematic pictures of the density profiles of the shell in (a) deceleration-dominated phase (Ss < S(joc)i (b) intermediate phase 
(Sjec < Es < 2Stioc)i E'lid (c) self-gravity-dominated phase {2Stioc < ^s)- 
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t/to 



Fig. 3. — Time evolution of the ratio of the colum n d ensity Eg 
to the characteristic column density S^cc (Equation l|16|l ^. Evolu- 
tionary phases corresponding to Figure [5] are labeled. 

sity Eg is given by the thin-shell approximation (see Sec- 
tion [2]). The column density Eg — pe-Rs/S can be de- 
rived from the mass conservation, since Mg = AttGR^'Ss 
(see Section [2]) . The time evolution of the density pro- 
file is roughly divided into the following three phases: 
deceleration-dominated phase (Eg < Edoc)j intermediate 
phase (Edoc < Eg < 2Edoc), and self-gravity-dominated 
phase (2Edcc < Sg), depending on the value of Eg/Edcc- 
The schematic pictures of the density profiles in these 
three phases are shown in Figures [2j Figure [3] shows 
the time evolution of Eg/Edec- In the early deceleration- 
dominated phase, Rc is outside of the shell and it is in 
front of the SF, or Rsf < Rc (see Figure d^a)). This 
means that the actual density peak exists at Rsf- As 
the shell expands, Eg increases by accretion while Edcc 
decreases by deceleration. In Figure [3l it is seen that 
Eg becomes larger than Edcc at t/to ~ 0.44. When 
Eg > Edcc, Rc is inside the shell. In the intermediate 
phase (Edcc < Eg < 2Edcc), Rc is closer to Rsf than 
i?CD as shown in Figure [IJb) . When Eg becomes larger 
than 2Edcc {t/to > 0.57), Rc becomes closer to i?cD than 
Rsf (see Figure [SJc)). Since the period of the inter- 
mediate phase is relatively short, roughly speaking, the 
density profile transforms from the deceleration- to the 
self-gravity-dominated profiles around t/to ^ 0.5. 



Obtained density profile by above semi-analytic 
method is compared with results of ID simulation. We 
use the ID sph erically symmetric Lagrangian Godunov 
method (jvan Le er 1997). We do not calculate the ra- 
diative transfer of ionizing photons and ionized gas, but 
the cold gas is pushed by interior pressure whose value 
is given by Equation The equation of state is as- 
sumed to be isothermal. We calculate the expanding 
shell around the 41M0 star that is embedded by the 
uniform ambient gas of uf — 10'^ cm~'^. 




-0.01 0.00 0.01 0.02 0.03 

(r-RcD)/Ro 

Fig. 4. — Snapshots of density profiles for t/to = 0.5, 0.7, 1.0, and 
1.26. The abscissae are the distance from the CD. The thick gray 
lines in the upper panel indicate the results of the ID simulation. 
The dashed lines in the upper panel represent the instantaneous 
hydrostatic density profiles. 

Figure |4] shows the snapshots of density profiles for 
t/to — 0.5, 0.7, 1.0, and 1.26. The thick gray hues repre- 
sent the results of the ID calculation. The dashed lines 
show the density profiles obtained from the semi-analytic 
method. It is seen that the semi-analytic method de- 
scribes the density profile of the ID simulation rea- 
sonably well. The density profile in the semi-analytic 
method is slightly lower than the results of the ID cal- 
culation because the actual CD expands a little slower 
than the mean radius of the shell Rs in the thin-shell ap- 
proximation. As shown in Section 13.11 one can see that 
the density peak moves the CD from the SF owing to the 
self-gravity. 



3.2. Comparison with One- Dimensional Simulation 



3.3. Scaling Law of Unperturbed Density Profiles 
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As shown in Section [21 the non-dimensional position of 
the sheU, Rg , is approximately independent of the model 
parameters (tie, Qvv)- Similarly, it is useful to inves- 
tigate how the density profile depends on the above pa- 
rameters. The non-dimensional pressures at the CD and 
the SF are given by Equation © and V^, respectively, 
that is, they are independent of the parameters. More- 
over, the pressure at the density peak Poo = PooCg is 
also independent of the parameters as seen in Equation 
(lisp . Thus, noting that the non-dimensional sound speed 
Cs — Csto/Ro is proportional to the reciprocal of the typ- 
ical Mach number Aio = 4i?o/(7toCs), where the factor 
of 4/7 arises from Equation (|4]), we have the scaling laws 

of Ho, Poo, and iff given by 



Ho cx clPoo^/' a Mo^ 
poo oc CjT^Poo oc Ml, 



and 



r —1/2 - n 
tff (x Poo « CsP[ 



-1/2 



00 



cx M 



' 



(20) 
(21) 



(22) 



respectively, where tff = l/V^vrGpoo is the free fall 
timescale of the shell. As a result, it is found that the 
density profiles for various set of (tie, Quv) are charac- 
terized by a single parameter, that is the typical Mach 
number, 



A^o = |^ = 7-^/^ 

7 Csto 

where T^^io = Tc/lO K. 



rp-l/2 -1/14 

UV,49-'c40 "e.S ' 



(23) 



4. INFLUENCE OF ASYMMETRIC DENSITY 
PROFILE ON GRAVITATIONAL INSTABILITY 

As shown in Section [3l the expanding shell has the 
highly asymmetric density profile and it is expected to 
influence the GI. In this section, we investigate influences 
of the asymmetric density profile on the dispersion rela- 
tion of the GI. What we discuss here is to extend the clas- 
sical stability analysis of the GI in the sym metr ic layer 
with respect to the mid-plane (iGoldreich fc Lvnden-Bell 



1965 : lElmegreen fc Elmegreen 1 119781 : iLubow fc Pringle 



1993[) to the GI in the asymmetric layer. The linear anal- 



ysis in the i ncompressible limit has been investigated by 
IVoit1(IT988l) . 

We take z-axis parallel to the thickness of the layer, 
and take a:-axis as the transverse direction. The density 
is assumed to peak at z = 0, and positions of bound- 
aries are zi and Z2 (zi < Z2). We consider a layer that 
is subject to a constant deceleration. The deceleration 
arises from the difference of pressures on the boundaries 
(z = zi and Z2). In this case, the position of the den- 
sity peak is not in the mid-plane of the layer, the density 
profile is asymmetric, or — zi ^ Z2. The amount of de- 
celeration directly enhances the degree of asymmetry of 
the density profile. 

4.1. Perturbation Equations 
We consider the following perturbations: 

piz,x,t)^poiz)+Spiz)e'^''--^'\ 
vAz,x,t) = v,{zy^'''--^'\ (24) 
v^{z,x,t) = v^{z)e^^^--'\ 



0(z, x, = 0o(^) + 5</)(z)e^('="-"*) . 
Perturbation equations are 

. . , d(po'U2) , ., „ 
— lujop H ^ h Poikvx — 0, (25) 



dz 



dz V Po 



ujVr. = A; I — 

PO 



and 



dz2 



k^5(j) = AnGSp, 



(26) 
(27) 

(28) 



where the sound speed is assumed to be constant. 

4.1.1. Boundary Conditions 

To concentrate on the effect of asymmetry of the un- 
perturbed state, we i mpose the CD boundary condi- 
tions at both zi and Z9 (IGoldreich fc Lynden-Bell II1965I : 
lElmegreen &: Eilmegreen II 19781) as follows: 

Sp{zi) = Szi, v^{zi) = -iujSzi, 

UZ 2 = Zi 



(W0 

dz 



and 



Sp{z2) = - 



dpo 



kS(f>{zi) + ATrGpo{zi)Szi, 



5z2, Vz{z2) — —iioSzi, 



dz Z = Z2 



+ kS(f>{z2) + 4:nGpo{z2)Sz2, 



(29) 
(30) 

(31) 
(32) 



where Szf and Sz2 are the displacements of the surfaces 
at zi and Z2, respectively. 

4.1.2. Numerical Methods 

We solve Equations (l^5t - (|28l) as a boundary-value 
problem for a given wavenumber. Equations (|25p . ()26|) . 
and (j28p are integrated from z = zi to Z2 by using the 
fourth order Runge-Kutta method. Note that Vx is de- 
termined by 5p and Scf) from Equation (j27p . Given ui, 
at z = zi, we have five unknown variables {5p, Vz, S(j), 
dS(t>/dz, and 5zi), and have three boundary conditions 
(see Equations ([59]) and pO|) ). Therefore, if we deter- 
mine two variables Qi and Q2, all variables at zi are 
specified, where Qi is one of {dp, Vz, Szi) and Q2 is one 
of {S(j>, dScfi/dz). Generally, the boundary conditions at 
Z2 are not satisfied if we start from arbitrary values of Qi 
and Q2 at Zf. Equation (|3ip can always be satisfied by 
using a linear combination of two independent solutions 
having the boundary values (Qi(zi), (52(^1)) = (l^O) 
and (0,1). Eigenvalue, uj, is modified iteratively until 
the solutions satisfy Equation (|32|) by using the Newton- 
Raphson method. 
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Fig. 5. — Dispersion relations for (a) the even mode and (b) the odd mode in the symmetric layer with \zi \ = 22 = 3Ho. The ordinate 
denotes ui^ /2nGpoo- The abscissa denotes the wavenumber multiplied by the effective thickness of the shell, H^g = cr/(2poo)- The dashed 
lines show higher harmonics of the sound wave with respect to the z-direction. 



O.i 




0.0008 



(b) SG-mode 



Fig. 6. — Distribution of the Lagrangian density perturbation, Ap/po = <5p/po + v^dln po{z)dz for kH^g = 3.5 in the even mode. Each 
panel corresponds to (a) the P-mode and (b) the SG-mode. The normalization is determined by \5z'2\/z2 = 0.1. 



4.2. Symmetric Layer 

First, we investigate the symmetric case with —Zi = 
Z2 — 3Hq. Because of symmetry, perturbation can be 
divided by even and odd modes completely. In the even 
(odd) mode, the density perturbation is symmetric (an- 
tisymmetric) with respect to the mid-plane. Figures 
[S]Ja) and (b) show the dispersion relations for the even 
and odd modes, respectively. The abscissa denotes the 
wavenumber multiplied by the effective thickness of the 
shell, Hcf[ = cr/(2poo)- It is well known that the unsta- 
ble mode (o;^ < 0) is found only in the even mode. The 
unstable mode belongs to the "compressible mode" that 
means that the density perturbation in the central re- 
gion collapses leaving behind the gas around boundaries. 
The detailed structure of stable modes is also plotted in 
Figure [H] The stable mode can be divided by the "P 
mode" (pressure mode, or compressible mode) and the 
"SG mode" (surface-gravity mode) . Figure [H] shows that 



the distribution of the Lagrangian density perturbation 
A/5 = Sp + Vzdpo/dz for kHcS = 3.5 in the even mode. 
Figures [Hfa) and (b) correspond to the P and the SG 
modes, respectively. The normalization is determined by 
\Sz2\/z2 = 0.1. One can see that Ap profiles of the P 
and SG modes are quite different. In the P mode, Ap/po 
peaks at the mid-plane. The displacement of the bound- 
ary \6z2\/z2 is negligible compared with Ap/po. The P 
mode propagates as longitudinal variation of pressure. 
On the other hand, the SG-mode has two Ap/po peaks 
near both boundaries, and has the minimum value at the 
mid-plane. Moreover, since \Sz2\/z2 is much larger than 
Sp/po, the SG mode is almost incompressible. The SG 
mode propagates as the deformation of the surface. In 
Figure ini^a), it is seen that the unstable mode transforms 
into the stable P mode around kHcs ~ 1. On the other 
hand, there is the another stable mode labelled by the 
SG mode {kHo < 2). The P and SG modes approach 
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each other as the wavenumber rises from the smaU Hmit. 
One can see a remarkable feature around fciJeff ^ 2.5 
where these two modes do not intersect but begin to 
move apart. At this point, they exchange their proper- 
ties, suggesting the mode exchange. The mode exchange 
between the P and the SG modes is also occurred in the 
odd mode (see Figure [5lb)). This is the first time when 
the mode exchange is found in the dispersion relation of 
self-gravitating layer. The dashed lines in Figures [5] show 
higher harmonics of the sound wave with respect to the 
z-direction. For large wavenumber, the frequencies of 
the P and the SG modes show different dependence on 
wavenumber. The P mode shows dependence while 
the SG mode shows k dependence. The angular frequen- 
cies of the SG modes associated by the deformation of zi 
and Z2 are 



and 



SG,Z2 



-2'kGpq{z2) 



1 dP^ 



Po dr 



1 dPn 



(33) 



(34) 



Po dr 

respectively (jWelter fc Schmid-Burgk 1119811) . for kHg > 
1. The SG branches for large wavenumber in the even 
and odd modes are identical with each other because the 
surface gravities at the two boundaries are the same. 

4.3. Asymmetric Layer 

In this section, we investigate the dependence of the 
dispersion relation on the degree of the asymmetry by 
changing zi . Since the layer is no longer symmetric with 
respect to z = 0, perturbations cannot be divided into 
the even and the odd modes. Figure [71[a) shows the dis- 
persion relation for zi = —2Hq. One can see more com- 
plex structure of the mode exchanges around kH^s ~ 2.4 
than that in Figure [S] For large wavenumber, the an- 
gular frequencies of the two SG modes split because 
'^SG zi < "^10 Z2- Figure IHKa) shows the cross section 
of the layer (zi = —2Ho) in the fastest growing mode. 
The contour indicates the density perturbation normal- 
ized by Poo- Here, we take (5pmax/poo = 0-2 to specify the 
normalization of the perturbations. The arrows repre- 
sent the velocity vectors. The boundary surfaces hardly 
deform, and the gas collapses from all directions to the 
center (z = 0, a; = 0). This behavior corresponds to the 
compressible mode. In Figure [71[a) , the unstable branch 
transforms the P mode around kH^s ^ 1 and it is con- 
nected with the SG mode around kHcs ^ 2.4 through 
the mode exchange. The case with stronger asymmetry 
with zi = —Hq is shown in Figure [TKb). In this case, 
the difference of the angular frequencies between the two 
SG modes is larger because the surface gravity at zi is 
lower. The frequency of the SG mode associated with 
zi becomes smaller than that with Z2- Comparing with 
Figure [TJa), as well as frequency, the wavenumber of the 
mode exchange is smaller. As a result, the frequency 
range of P mode is narrower and the range of the SG 
mode spreads. The P mode oc k^ is expected to dis- 
appear when the wavenumber of the mode exchange is 
smaller than a critical wavenumber that separates unsta- 
ble mode from stable mode. 

The case with —zi< Hq is quite different from the case 
with —zi > Hq. The dispersion relation for —zi = O.SHq 



is shown in Figure [TKc). In Figure [7Kc), one can see that 
the angular frequency of the SG mode w|q ^.^ is signifi- 
cantly lower than WgQ for large wavenumbers. Unlike 
the case with —zi> Hq, the unstable mode appears to 
directly connect with the SG mode around kH^s 1.2 
as mentioned above. Figure [5fb) shows the cross section 
of the layer (zi = —O.SHq) in the fastest growing mode. 
The gas tends to collect toward the density peak z = 
because the unperturbed gravitational potential has the 
minimum value there. One can see that eigen-functions 
in z > are similar to those in Figure [D^a). The gas col- 
lapses toward the center (z = 0, a; = 0) leaving behind 
the gas around Z2. However, in the region where z < 0, 
eigenfunctions are quite different. The sound wave can 
travel between zi and the density peak many times dur- 
ing the development of the GI. Therefore, collapse toward 
z = is suppressed by the pressure gradient. However, 
the GI can proceed even in z < through the defor- 
mation of the zi that makes the gravitational potential 
deeper. From Figure |8l[b) , one can see that the veloc- 
ity field is not headed for the density peak (z = 0) but 
arises so that the surface at zi deforms. Therefore, the 
features of GI in the region z > and z < have proper- 
ties of "compressible mode" and "incompressible mode" , 
respectively. If the distance of the zi from the density 
peak is zero, the layer is unstable for all wavenumbers 
(see Figure EKd)). 

Figure [9] shows the dispersion relation of the unstable 
mode for variety of zi with Z2 = SHq. The thick solid 
and the thick dashed lines correspond to —zi/Hq = 2 
and 1, respectively. For —zi/Hq > 1, the growth rate 
— a;^/27rGpoo decreases as —zi/Hq decreases. This prop- 
erty is the same as that in sy mmetric layers (e.g., see 
Figure 1 in lNagai et al. IITQQSI ). On the other hand, for 
the cases with —zi/Hq = 0.3 and 0.1, Figure [9] shows 
that the maximum growth rate increases as —zi/Hq de- 
creases, indicating the opposite tendency to the case with 
—Zi/Hq > 1. On the other hand, the wavenumber of the 
most unstable mode is not different so much. One can 
see that the square growth rate in large wavenumber is 
proportional to cx fc while the square growth rate for 
—Zi/Hq > 1 is proportional to oc fc^. The unstable mode 
appears to directly connect with the surface gravity mode 
at zi whose frequency is given by Equation (|33p . This is 
also seen in Figure fTOl for k > fcmax- For —zi /Hq = 0, 
the surface gravity mode at zi becomes unstable for all 
wavenumber. In this case, destabilized surface gravity 
wave has the growth rate — — 27rGp(zi) < in- 

dependent of k for large wavenumber limit (see Equa- 
tion ((33)l ). This is the c ase with a static she l l with 
gdoc = (Equation iTomisaka fc Ikeuchi I (|1983D 

investigated this situation including shell curvature and 
fou nd th at the sh ell is unstable fo r all wavenumber (also 
see l Welter fc' Schmid-Burgk llT98lh . For-zi/i7o = -0.1, 
the square growth rate increases as oc fc with wavenum- 



bers because 



1 dPg 

Po dr 



fc < (Equation (|TT|) V 

2 

This is well-known scaling law of the Ray leigh- Taylor 
instability. The enhancement of the growth rate for 
—Zi/Hq < 1 arises from the combination of the GI and 
the Rayleigh- Taylor instability. 




Fig. 7. — Dispersion relations for —z\/Ho =(a)2.0, (b)l.O, (c)0.3, and (d)O.O. The ordinate and the abscissa are the same as Figure|5] 



5. GRAVITATIONAL INSTABILITY OF 
EXPANDING SHELLS 

In previous section, we focus on the effect of asymme- 
try of the density profile by imposing the same boundary 
conditions in both boundaries. In this section, in a more 
realistic situation, we investigate the stability of expand- 
ing shells driven by the expansion the HII region. The 
unperturbed density profile at each instant of time is 
given by the semi-analytic method presented in Section 
[21 We neglect the curvature effect and solve the per- 
turbation Equations (l25t - (|28l) but z — >■ r. In this case, 
as well as in the asymmetric density profile, the differ- 
ence of boundary properties between leading (the SF) 
and trailing (the CD) surfaces plays important roles in 
the GI. 

5.1. Influences of Boundaries on the Gravitational 
Instability 

Before presenting a linear analysis, we review how the 
SF and the CD influence the GI through the boundary 
effect. This point is important in understanding the re- 
sults of the linear analysis. In the early phase, the shell 
is highly confined by the ram pressure on the leading 
surface and by the thermal pressure on the trailing sur- 
face. In this phase, the pressure at boundaries is as 
large as that at the density maximum, and the thick- 
ness of the shell is much smaller than the scale height. 



-ffo = Cs/a/SttGpoo, where poo is the maximum density. 
Thus, in this phase, the boundary effect can strongly 
infiuence the GI. In the later phase, the boundary ef- 
fect of the CD is expected to be also important because 
the density peak is close to the CD as shown in Section 
13.11 In this section, we summarize how the growth rate 
of GI is controlled by the different boundary conditions. 
For simplicity, in Section 15. 1[ the layer is assumed to be 
symmetric with respect to the mid-plane, and physical 
variables are averaged across the thickness. 

5.1.1. Shock- confined Layer 

Many authors have investigated influences of the SF on 
the GI (Vishniac 1983: Elmegreen . igS^ iMshTI [l992l: 



Vishniac I 119941: 



19941: iWhitworth et all 



. - - ■ .Elmegreen I . 

1994at llwasaki fc Tsuribe 1 12008| ). The dispersion rela- 
tion of the shock-confined layer is given by 



2 I 2 



(35) 



where Eg is the column density, and k is the trans- 
verse wavenumber of the perturbation. This disper- 
sion relation is the same as that for the infinitesi- 
mally thin layer. In the highly confined layer, it 
is well known that the perturbation behaves like in- 
compressible mode because the sound-crossing time 
over the thickne ss is much smaller than the free-fall 
timescale, ~ 1/VGpoo (,Elmegrecn fc Elmegreen 1 119781 : 
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Fig. 8. — Cross sections of the layers with the fastest growing 
mode for (a)— 21 = 2.0//o and (b) — zi = 0.3-ffo. The contour 
indicates the density perturbation normalized by poo . The contour 
levels of the density perturbation take values of 0.04, 0.08, and 
0.16. 
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Fig. 9. — Dispersion relation of the asymmetric layer for 
—z\/Hq = 2(the thick solid line), l(the thick dashed line), 0.3(the 
thick dotted line), 0.1 (the thin solid line), 0(the thin dashed line), 
and —0.1 (the thin dotted line), where Z2/H0 is assumed to be 3.0. 



ILubow fc Pringle1ll993[) . Therefore, density fluctuation 
is small. The layer becomes unstable mainly by the de- 
formation of the surfaces that makes the perturbation 
of the gravitational potential deeper. Hereafter, we call 
this mode the "incompressible mode" . The deforma- 
tion of the SF generates the tangential flow carrying the 
gas from the convex to the concave regions (seen from 



the downstream). Therefore, the tangential flow tends 
to make the SF flat, suggesting that it suppresses the 
growth of the GI. In the shock-confined layer, the restor- 
ing term c^k'^ arises from the tangential flow behind the 
oblique SF. On the other hand, in the case of the in- 
finitesimally thin layer, this term comes from the pres- 
sure gradient. Therefore, the origin of the restoring force 
is quite different. From Equation (I35p . the maximum 
growth rate is given by ttGSs/cs, and the correspond- 
ing wavenumber is given by TrGSg/Cg. When is small 
{"^ PoqHq), the maximum growth rate is smaller than the 
inverse of the free-fall timescale, i/Gpoo, and the corre- 
sponding scale is larger than the scale height Hq that is 
comparable to the Jeans scale. 

5.1.2. Pressure- confined Layer 

Next, we review the influence of the CD on the GI. We 
consider the layer confined by thermal pressure of hot 
rarefied gases (CD boundary condition) on both sides. 
The dispersion relation becomes 

uj^ ~ 2nGJ:sLsk^ 



27rGfcEs 



(36) 



where Lg is the thickness of the layer, and we con- 
sider the large-scale limit where k <C 1/is- Detailed 
de rivation o f Equation (1361) is shown in Appendix 3 
of llwasaki fc Tsuribe . (2008.). In the pressure-confined 
layer, the stabilization effect of the tangential fiow does 
not exist. Therefore, the restoring term in Equation p6l) 
is quite different from that in Equation ([55)1 . Using the 
gravitational acceleration at the surfaces as \g\ — 27tG'Ss, 
we can express the restoring term in Equation (j36p by 
\g\Lsk'^. Therefore, one can see that the restoring force 
arises from the surface gravity wave. The layer with the 
CD boundary condition is less stabilized compared with 
the shock boundary condition because the phase veloc- 
ity of the gravity wave y'Jg\Ls is much smaller than Cg 
when Ls ^ H^. From the dispersion relation ([55]) . the 
maximum growth rate is comparable to the inverse of 
the free fall time of the layer ~ y'Gpoo, and the corre- 
sponding wavelength is about the thickness of the layer, 
c± L s (Elmcgrecn & Elmcgrccn 1978; Lubow & Pringle] 
|1993( ). Therefore, one can see that the most unstable 
mode in the pressure-confined layer has a larger growth 
rate and a sma ller scale than the shock -confined layer 
(see Figure 9 of llwasaki fc Tsuribe II2008L in detail). 

5.1.3. Expanding Shells 

Equations (|35l) and (p6| cannot be applied directly in 
the GI of the expanding shells because the GI is expected 
to be stabilized by evolutionary effects, such as the ex- 
pansion of_yieshe]landthe accretion of fresh gas through 
the SF. lElmegreen I ()1994[ ) derived the following approx- 
imate dispersion relation. 



(37) 




21,2 



The terms with Vg/Rs come from evolutionary effects 
that stabilize the GI. 

One can see that Equation ([57)) for the limit of 
Vs/Rs — > is th e sa me as Equation (1551). Th erefore, 
lElmegreenI (|1994[ ) and 'Whit worth et al. I |l994bD essen- 
tially applied Equation ([35)) in the context of the GI of 
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the expanding shell. However, they did not take into 
account the boundary effect of the CD on the trailing 
surface. Comparing Equations (|35)) and (p6)) , we suggest 
that the stability of the thin shell neglecting the effect of 
the CD is suffered by large stabilizing effect, and it will 
underestimate the growth rate of GI in expanding shells. 

5.2. Boundary Condition 

First, we assume that a constant pressure ex- 
erts on the CD all the time (the CD b ound- 
ary condition; iGoldreich fc Lvnden-Bell I 119651 : 
lElmegreen fc Elmegreen I I1978D . The boundary condi- 
tions are 

„^ 
dr r=flcD 

(38) 

and 

d5(j) 
dr 



5p{R 



CD J 



(5i?CD, Vr{RcT)) = ibjSRci 



kS4> + ATrGp{RcB)SRci) = 0, 



(39) 



where SRcn is the displacement of the CD. 

Next, let us consider the boundary conditions at r = 
i?SF- Since the unperturbed state is assumed to be 
the hydrostatic configuration, it is impossible to impose 
the shock boundary conditions self-consistcntly. In or- 
der to treat it self-consistently, time-dependent initial 
value problem is needed to be solved (jWelter 1 119821 : 
llwasaki fc Tsuribe I [20081) . Therefore, in this paper, we 
mimic the shock boundary conditions by introducing the 
stabilization effect. We consider the following two ap- 
proximate boundary conditions. 

Rigid surface boundary condition (RSBC). iVoit I 

(11988) and lUsami. Hanawa fc Fuiimoto I (|1995[ ) assumed 
that no ripples arise on the surface, or i5i?sF = 0, where 
(5i?SF is the displacement of the SF. The reason why we 
adopt (5i?sF = is that the thin-shell linear analysis of 
the layer confined by rigid surfaces gives the same disper- 
sion relation as that of the shock-confined layer (Equa- 
tion p5p ). In more precisely, in the shock-confined layer, 
the tangential flow boosts the suppression effect against 
the self-gravity as mentioned in Section 15.1.11 Instead, 
the RSBC weakens the self-gravity. 

Tangential flow boundary condition (TSBC). — If the 
SF is rippled, the tangential flow behind the SF is gen- 
erated. Therefore, we set the tangential velocity Vx at 
i?SF- Linearizing the Rankine-Hugoniot relation, we 



r 

have 



Vx{Rsf) = - ^SF - 



R 



SF 



ik5R. 



SF- 



(40) 



The detailed deriv a tion of Equation pO|) is found in 
llwasaki fc Tsuribe I (|2008l ). 

With both of above boundary conditions (RSBC and 
TFBC), we also impose the following ordinary used 
boundary conditions, 

Wr(i?SF) = iOjSRiiY, (41) 



and 



dr 



k5(j) + ATTGpiRsF)5RsF = 0. (42) 



It is well known that the SF of the deceler ation shell 
is subject to the hydrodynamical overstability (jVishniac I 
[T983i) . The linear analysis in this paper cannot capture 



the Vishniac instability (VI) correctly since the approxi- 
mate shock boundary conditions are imposed. The effect 
of the VI is discussed in Section |6l 
The numerical method is the same as that in Section 

5.3. Scaling Law of Dispersion Relations 

As shown in Section 13.31 it is found that the den- 
sity profiles are characterized by a single parameter Mq. 
This is because the scale height, the peak density, and 
the free fall time have the scaling laws with respect to 
A^o as shown in Equations (|20)) - (|22|) . The same is the 
case with the perturbation equations and the disper- 
sion relation. The non-dimensional maximum growth 
rate Wmax = i^max^o and the corresponding wavenumber 
fcmax = fcmax-Ro scalc as oc ^Ao and cx A4q, respectively. 
Therefore, in the present model, the evolution of the shell 
for various set of (tt-e, Quv) can be described by a single 
unperturbed profile and a single time-dependent disper- 
sion relation that are normalized by Hq , poo i and tg . The 
result can be applicable to a wide range of parameters 
simply by using the scaling relation on Mq. 

5.4. Results 




Fig. 10. — Dispersion relations derived from our linear analysis 
using RSBC (the solid lines) and TFBC (the dotted lines). For 
com paris on, the dispersion relation of shock-confined layer (Equa- 
tion l|35|l ) is plotted by the dashed gray lines. The thick gray lines 
represent modified dispersion relation (Equation Il43l l). The ab- 
scissa and ordinate axes indicate the wavenumber kHo/2n and the 
growth rate normalized by iff = 1/\/2ttGpoo, respectively. 



At any time, the unperturbed state is given by the 
procedure in Section |31 Perturbation Equations (|25p - 
(128]) are solved as the eigenvalue- and boundary-value 
problem. As a result, the growth rate, uj{k,t) can be 
obtained as a function of the wavenumber and time. 

First, we present the results of the linear analysis in 
Figure [10] at various epochs. The ordinate and the ab- 
scissa axes represent the non-dimensional growth rate 
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Fig. 11. — Time evolution of the ratio of Cgg to Cs. 
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Fig. 12. — Cross section of the shell is plotted using eigen- 
functions. The corresponding time and angular wavenumber are 
t/to = 1.3 and Z = 52, respectively. The contour indicates the den- 
sity perturbation normalized by poo • The vectors represent velocity 
perturbations. 

ujtg and wavenumber fc77o/27r. The solid and the dot- 
ted Unes indicate the results of the linear analysis us- 
ing RSBC and TFBC, respectively. We refer the growth 
rates obtained by using RSBC and TFBC to wrsbc and 
wtfbc? respectively. The dependence of the dispersion 
relation on the parameters (tie, Quv) can be eliminated 
by using non-dimensional growth rate Lotg and wavenum- 
ber kHo as shown in Section 15.31 We have confirmed 
that the dispersion relation is identical to that with 
other parameter sets of (tie, Quv) by using the non- 
dimensional quantities. Figure ITOl shows that the differ- 
ence between wrsbc and wtfbc is negligible although 
RSBC and TFBC are physically quite different. 

In this analysis, we do not take into account the evo- 
lutionary effects, such as the expansion and accretion of 
the gas. Therefore, we compare the results of the lin- 
ear analysis with the dispersion relation of the shock- 
confined layer (Equation (|35|)) rather than that of the 
expanding shell (Equation (157)) ). One can see that the 
growth rate is larger than the prediction from the shock- 
confined layer. As shown in Section \57\\ this difference 
comes from the boundary effect of the CD. Therefore, 
the shell is expected to begin to grow earlier and to frag- 
ment more quickly than the prediction from lElmegreenl 
(|1994[) that is based on Equation ([55]) . 

The dispersion relation with CD + SF boundary con- 
ditions is expected to lie between that with SF -I- SF 
(Equation ^) and that with CD + CD (Equation 



(|36|) V Therefore, to approximate the dispersion relation 
with RSBC analytically, we combine Equation (|35|) with 
Equation (|36l) as follows: 

(43) 

where CcS is the effective sound speed, 



Ceff = A27:GJ:,L,s +(!)', (44) 

where A is a parameter and L^s ~ Es/poo is the effec- 
tive thickness that approaches the actual thickness Ls 
for small Sg and 2Ho for large Sg- The first and the 
second terms inside the square root correspond to the ef- 
fect of the CD and the SF boundary conditions, respec- 
tively. Here, we choose the parameter A by the condition 
where the maximum value of Wmod coincides with that 
of wrsbc- As a result, it is found that a single value of 
A = 0.39 shows good agreement in growth rates between 
the modified dispersion relation and the detailed linear 
analysis. The modified dispersion relations in Equation 
([43|l are plotted by the thick gray lines in Figure [101 In 
Figure [TOl one can see that Wmod well describes wrsbc 
for k < fcmax all the time. This suggests that simply 
Equation can describe the most unstable mode ob- 
tained by the detailed linear analysis all the time. Figure 
[11] shows the time evolution of the effective sound speed. 
In the early phase, CcS is about 0.5cs, suggesting that 
the effect of the CD diminishes the effective sound speed 
Ceff by half in Equation (j35|) . As the shell expands, c^s 
increases. 

The effect of asymmetric density is seen in Figure [TUj 
where it is found that Wmod deviates from wrsbc for 
k > fcmax in the gravity-dominated phase (t/to > 0.5). 
This is because oj'^^^ connects with the P mode cx fc^ 
while t^^sRC connects with the SG mode cx fc as shown 
in SectiongJ 

The predicted cross section of the shell by the linear 
analysis for (Quv = lO"®-^^ s~\ tie = 10^ cm"^) is 
shown in Figure [T^ by using the eigenfunctions. The cor- 
responding time is t/to = 1.3 and the angular wavenum- 
ber is ^ = 52. In Figure [T^ the gas tends to accumulate 
onto the peak only through the upper half region r > R^. 
This property of the fiow can be seen from the direction 
of arrows in Figure 1121 Actually, in the upper half re- 
gion, we find i?sF — Rc = 1.05i?o > Hq that represents 
that the gas can collapse to the peak because the sound 
wave cannot travel from i?c to i?sF within the free fall 
time. On the other hand, in the region of bottom half 
[r < Rc), we find that i?c - Rcu = 0.2857fo < Ho- This 
indicates that the gas in r < i?c cannot collapse to the 
peak because the sound wave can travel from Rc to Rcd 
many times within the free fall time. Thus, the pressure 
gradient prevents the compression of gas in the region 
r < Rc- However, the GI can proceed even in r < i?c 
through the deformation of the CD that makes the grav- 
itational potential deeper. Therefore, the features of GI 
in the region r > Rc and r < Rc have the properties 
of the "compressible mode" and "incompressible mode" , 
respectively. 

6. DISCUSSION 

The gravitational fragmentation of expanding shells 
confined from both sides by the CD was investigated 
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bvlDale et al.l ()2009D numerically and bv lWiinsch et al. 1 
(|2010D using analytical approximations. They assumed 
that the thermal pressure on both sides is the same and 
temporally constant. Therefore, the density peak is al- 
ways around the mid-plane of the shell, and the den- 
sity profile is almost symmetric. In their calculation, 
the column density decreases with time because the shell 
expands keeping the mass fixed. Therefore, the pres- 
sures at the boundaries approach to the peak pressure. 
They found that the confining pressure accelerates frag- 
mentation in the later phase, and described this effect 
as "pressure-assisted" gravitational fragmentation. This 
mod e is the same as t he inc ompressible mode in this pa- 
per. iWiinsch et al. I (|2010D establish ed a semi-analyti c 
linear analysis that explains results of lDale et al.l ()2009[ ). 

The linear analysis in this paper cannot describe the 
VI correctly since the approximate shock boundary con- 
ditions are impo sed in Section[5l The original analysis by 
iVishniad ()1983l) did not find the finite scale most unsta- 
ble mode because the thickness of the shell is neglected. 
iVishniac fc Rvul (|l989) derived a simple analytic dis- 
persion relation of the VI for a decelerating isothermal 
spherical shock wav e taking into account th e effect of the 
thickness (also see iRvu fc Vishniac I Il987| ) . Although, 
their analysis did not include the self-gravity, here, we 
use their dispersion relation (see Equations 19(a) and 
(b) in their paper) to estimate the effect of the VI. Their 
dispersion relation depends on the Mach number M of 
the shell and the expansion law. For the case with the 
expanding HII regions, the shell expands as oc t"*/^ if the 
self-gravity is neglected. In this case, the perturbation 
grows not exponentially but in a power-law cx t", where 
s characterizes the growth rate. Figure [T51 shows the real 
part of s as a function of the angular wavenumber I. One 
can see that the maximum growth rate Re(s) increases 
with Ai ■ The angular scale of the most unstable mode 
is smaller for larger A4 . We find that the unstable mode 
exists only for J\4 > 4.7. To see the typical value of the 
Mach number, we consider the expanding shell around 
the 4IM0 star that is embedded by the uniform ambient 
gas of riE — 10'^ cm^'^. FigurefMlshows the Mach number 
of the sheU for Tc = 10 K (the sohd line) and 30 K (the 
dashed line) . In the early phase when the self-gravity is 
not important {t/to < 0.5), since the Mach number is 
as large as several tens, Re(s) is large. The small scale 
perturbation with I = 10^ ^ 10^ quickly grows and satu - 
rates in the nonlinear stage (|Mac Low fc Norman |[l993l ) . 
On the other hand, in the later phase (the self-gravity- 
dominated phase, t/to > 0.5), the Mach number is as low 
as 5 — 10 as shown in Figure [T4l In this phase, Re(s) ~ 1 
from Figure 1131 This means that the growth rate of the 
perturbations is comparable to the expansion rate cx t^/*". 
Therefore, in the self-gravity-dominated phase, the VI is 
not expected to be important. The infiuence of VI on 
the GI is expected to be only the increase of the initial 
amplitude of perturbations for the GI. 

7. SUMMARY 

In this paper, we have performed linear perturbation 
analysis of decelerating shells created by the expansion 
of HII regions. We summarize our results as follows: 

1. We develop a semi-analytic method for describing 
the density profile in the shell. The time evolution 
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Fig. 13. — Growth rate of the VI when the shell expands as 
oc t^/^. Each line corresponds to Ai = 4.7, 10, 7, and 4.7. 




Fig. 14.— The Mach number of the shell for Quv = 10"'^ s"^ 
and nE = 10'^ cm~^. The solid and the dashed lines indicate the 
case with Tc=10 K and 30 K, respectively. 

of the density profile of the expanding shell can be 
divided into three phases, deceleration-dominated, 
intermediate, and self-gravity-dominated phase. In 
the deceleration-dominated phase, the density peak 
is in SF by the inertia force owing to the deceler- 
ation. As the shell mass increases and the self- 
gravity becomes important, the density peak is in- 
side the shell, but it is closer to the SF than the 
CD in the intermediate phase. In the self-gravity- 
dominated phase, the shell becomes massive and 
the density peak is closer to the CD than the SF. 
The evolution is confirmed by ID hydrodynamical 
simulation. 

2. We show detailed structures of dispersion relation 
in the asymmetric layer subjected to a constant 
deceleration both of unstable and stable modes by 
imposing the CD boundary condition from/at both 
sides. 

• We discover the mode exchange between the 
compressible and surface-gravity modes in the 
stable regime. 

• In a situation where the distance from one sur- 
face zi to the density peak z = is smaller 
than the scale height of the self-gravity Hq 
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and the distance from the other surface zi to 
z = larger than _ffoi the nature of the GI is 
quite different from the symmetric case with 
the same cohimn density and the peak density. 
The eigenfunction in the region < z < Z2 
is approximately the compressible mode. On 
the other hand, the eigenfunction in the re- 
gion zi < z < is approximately the incom- 
pressible mode. Moreover, the growth rate 
is enhanced compared with symmetric cases 
through cooperation with the Ray leigh- Taylor 
instability. 

3. We investigate linear stability of expanding shells 
driven by HII regions taking into account the 
shock-like boundary condition on the leading sur- 
face, the CD boundary condition on the trailing 
surface, and the asymmetric density profile ob- 
tained by the semi-analytic method. 

• The shell is expected to grow earlier than 
the prediction of previous studies ([Elmegreen I 



ITOOllWhitworth et al. l[T994bl ) that are based 
on the dispersion relation of the shock- 
confined layer. 

• In the self-gravity-dominated phase, since the 
density peak is closer to the CD than the SF, 
the CD is expected to deform significantly. 

These results provide useful knowledge for the analysis 
of more detailed nonlinear n umerical simulations that is 
the scope of our next paper (jlwasaki et al. |[2011|) . 

We thank the referee for many constructive comments 
that improve our paper significantly. This work was 
supported by Grants-in-Aid for Scientific Research from 
the MEXT of Japan (K.I.:22864006; S.I.:18540238 and 
16077202), and Research Fellowship from JSPS (K.I.:21- 
1979). Thi s work was based on the results of the compan- 
ion paper (jlwasaki et al. ir201lD where numerical compu- 
tations carried out on Cray XT4 at the CfCA of Numeri- 
cal Astronomical Observatory of Japan. The page charge 
of this paper is supported by CfCA 
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